#ifndef BS_BOS_RESIDUAL_2D
#define BS_BOS_RESIDUAL_2D

#include "boundary_building.h"
#include "fadiff.h"
#include "residual1D.h"
#include "coefficient.h"

namespace blue_sky{

	template<class strategy_t>
	class Residual2D
	{
	public:

		typedef typename strategy_t::csr_matrix_t           csr_matrix_t;
		typedef typename csr_matrix_t::fp_vector_type		fp_vector_type;
		typedef typename strategy_t::index_t				index_t;

		template <class T>
		T operator()(const T& U_ij, const T& U_ipj,const T& U_imj, const T& U_ijp, const T& U_ijm, 
			const int cell, double dt, smart_ptr<build_boundary<strategy_t>, true> boundary_builder)
		{
			Residual1D<strategy_t> residual1D;
			Coefficient<strategy_t> a;
			fp_vector_type area = boundary_builder->geom_factors;
			return residual1D(U_ij, U_ipj, U_imj, cell, dt, boundary_builder) 
				+ dt * ( 
				- area[cell * 4 + Y2] * a(U_ij) * (U_ijp - U_ij)
				+ area[cell * 4 + Y1] * a(U_ij) * (U_ij - U_ijm)
				);
		}
	};

	template <class C, class strategy_t>
	class ResidualDiff2D
	{
	public:

		template <class T>
		T operator()(const int cell, double dt, const T& U_ij_, const T& U_ipj_,const T& U_imj_, const T& U_ijp_, const T& U_ijm_, T& dRdUij, T& dRdUipj, T& dRdUimj, T& dRdUijp, T& dRdUijm)
		{
			// Initialize arguments
			F<T> U_ij(U_ij_), U_ipj(U_ipj_),
				U_imj(U_imj_),U_ijp(U_ijp_),
				U_ijm(U_ijm_); 

			//difirentiation
			U_ij.diff(0,5);
			U_ipj.diff(1,5);
			U_imj.diff(2,5);
			U_ijp.diff(3,5);
			U_ijm.diff(4,5);

			C func;             // Instantiate functor
			F<T> f(func(U_ij, U_ipj, U_imj, U_ijp, U_ijm, cell, dt, boundary_builder_));  // Evaluate function and record DAG

			//save derivatives
			dRdUij=f.d(0);
			dRdUipj=f.d(1);
			dRdUimj=f.d(2);
			dRdUijp=f.d(3);
			dRdUijm=f.d(4);

			return f.x();       // Return function value
		}

		ResidualDiff2D(smart_ptr<build_boundary<strategy_t>, true> boundary_builder){
			this->boundary_builder_ = boundary_builder;
		}

		smart_ptr<build_boundary<strategy_t>, true> boundary_builder_;
	};

}
#endif //RESIDUAL